library(survival)
geno<-read.table("input_expr",sep="\t",head=T)
ordersamples<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/ordered_sampes_Martingale.txt")
ordersamples<-as.vector(ordersamples[,1])
ordersamples<-gsub("-",".",ordersamples)
fil_snps<-read.table("/data1/bsi/BORA_processing/devel/eqtl/permu_geneexpr_topgenes/uniq_genes_maintop8000.txt",sep="\t")
fil_snps<-as.vector(t(fil_snps))
geno<-geno[as.character(fil_snps),]
i<-1
hazardrank<-matrix(nrow=nrow(geno), ncol=ncol(geno))
colgeno<-colnames(geno)
colnames(hazardrank)<-colgeno
clin<-read.delim("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/clinical" ,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="",
header=TRUE,
col.names=c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE"),as.is=c(TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),colClasses=c("character","numeric","factor","numeric","numeric","numeric","integer","integer","integer","factor"))
clin$SUBJECT<-gsub("-",".",clin$SUBJECT)
samplesize<-20
num<-20
while(i<nrow(geno)+1)
{
	geno1<-geno[i,]
	d<-1
	while(d<num+1)
	{
		s<-d;
		sub<-c()
		main<-c()
		if(d > 1)
		{
			main<-seq(1,d-1,1)
		}
		nu<-1
		while(s<length(geno1)+1)
		{
			sub[nu]<-s
			sequ<-c()
			if(s+1 <= length(geno1))
			{
			sequ<-seq(s+1,s+num-1,1)
			if(s+num-1 > length(geno1))
			{
			sequ<-seq(s+1,length(geno1),1)
			}
			}
			main<-c(main,sequ)
			nu<-nu+1
			s<-s+num
		}
		genomain<-geno1[main]
		genosub<-geno1[sub]
		dt <- merge(clin,t(genomain),by.y="row.names",by.x="SUBJECT")
		colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		j=1
		order_samp<-c()
		while(j<length(ordersamples))
		{
			order_samp<-c(order_samp,which(dt[,1]==ordersamples[j]))
			j<-j+1
		}
		order_samp<-c(order_samp,which(dt[,1]==ordersamples[length(ordersamples)]))
		dt<-dt[order_samp,]
		fit1 <- coxph(Surv(OS,ALIVEOS) ~ EXPR+AGE,data=dt)
		s<-summary(fit1)$coefficients
		e_signif<--log10(as.numeric(s["EXPR","Pr(>|z|)"]))
		effect<-as.numeric(s["EXPR","coef"])
		e_signif__effect<-paste(e_signif,effect,sep="__")
		dt1 <-merge(clin,t(geno1),by.y="row.names",by.x="SUBJECT")
		colnames(dt1) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		j=1
		order_samp<-c()
		while(j<length(ordersamples))
		{
			order_samp<-c(order_samp,which(dt1[,1]==ordersamples[j]))
			j<-j+1
		}
		order_samp<-c(order_samp,which(dt1[,1]==ordersamples[length(ordersamples)]))
		dt1<-dt1[order_samp,]
		pred <- predict(fit1, dt1,type="lp")
		wanted_ranks<-rank(pred)
		hazardrank[i,sub]<-wanted_ranks[sub]
		d<-d+1
	}
	print(i)
	i<-i+1
}
deltarank<-hazardrank
rownames(clin)<-clin[,1]
clin1<-clin[colnames(hazardrank),]
fit<-coxph(Surv(OS,ALIVEOS)~AGE,data=clin1)
pred<-predict(fit,clin1,type="lp")
d<-rank(pred)
surv<-survConcordance(Surv(OS,ALIVEOS)~AGE,data=clin1)
i<-1
while(i<nrow(hazardrank)+1)
{	
	deltarank[i,]<-hazardrank[i,]-d
	i<-i+1
}
rownames(hazardrank)<-rownames(geno)
rownames(deltarank)<-rownames(geno)
i<-1
deltaAUC<-c()
while(i<nrow(deltarank)+1)
{
	clin2<-cbind(clin1,deltarank[i,])
	colnames(clin2)[11]<-"diff"
	fit<-coxph(Surv(OS,ALIVEOS)~diff,data=clin2)
	surv<-survConcordance(Surv(OS,ALIVEOS)~diff,data=clin2)
	deltaAUC<-c(deltaAUC,surv$concordance)
	if((surv$concordance-1.96*surv$std.err)>0.5)
	{
		array<-c(array,rownames(hazardrank)[i])
		print(rownames(hazardrank)[i])
	}
	i<-i+1
}
deltaAUC<-cbind(rownames(geno),deltaAUC)
colnames(deltaAUC)<-c("marker","deltaAUC")
write.table(deltaAUC,"deltaAUC.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
pdf("heatmap_afterfilter.pdf")
heatmap(hazardrank[as.character(array),],cexRow=0.5,cexCol=0.25)
heatmap(deltarank[as.character(array),],cexRow=0.5,cexCol=0.25)
dev.off()
write.table(hazardrank,"hazard.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
